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NUMERICAL  SIMULATIONS  OF  ASYMMETRIC  MIXING 
IN  PLANAR  SHEAR  FLOWS 

I.  INTRODUCTION 

This  paper  describes  numerical  simulations  of  shear  flows  and  the 
development  and  evolution  of  ooherent  structures.  There  are  two  aspeots 
of  this  problem  which  are  addressed  here.  The  first  is  developing  the 
numerical  model  that  was  used  in  these  studies.  In  particular,  we  are 
concerned  with  the  treatment  of  inflow  and  outflow  boundary  conditions 
suitable  for  both  compressible  and  incompressible  flows.  The  second 
aspect  is  using  this  model  to  describe  shear  flows  in  a  splitter-plate 
configuration. 

The  problem  of  developing  the  proper  computational  tools  has  been 
discussed  previously  (Boris  et  al.  1985).  He  elaborate  on  this  subject 
in  Section  III.  The  numerical  model  we  use  is  a  restructured  version  of 
the  FAST2D  computer  code.  This  incorporates  the  Flux-Corrected  Transport 
(FCT)  continuity  equation  algorithm  (Boris  1976b,  Boris  &  Book  1976)  which 
has  been  tested  extensively  for  shook,  blast,  detonation,  fluid  instability 
and  beam-generated  turbulence  calculations  (e.g.,  Book  et  al.  1980,  Oran 
et  al.  1982,  Pioone  &  Boris  1983).  Since  the  algorithm  is  explioit,  the 
code  is  best  for  studying  flows  that  move  at  a  substantial  fraction  of  the 
speed  of  sound  in  the  material.  Using  time  step  splitting  techniques,  FCT 
may  be  coupled  to  algorithms  representing  other  physical  processes,  such 
as  diffusion  and  conduction  (Oran  and  Boris,  1981). 

The  second  aspect,  application  of  the  model,  is  the  main  goal  of  this 
work.  He  have  used  the  model  to  simulate  time-dependent  flows  in  the 

splitter-plate  configuration  for  whiob  substantial  data  exist  on  the 
Manuscript  approved  May  28, 1985. 


transition  to  turbulence.  Our  interest  here  is  in  the  entrainnent  and 
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mixing  in  the  resulting  ooherent  structures  in  the  flow.  In  Section  IV, 
we  describe  the  computational  results  and  their  implications. 


KKa 


II.  BACKGROUND 


The  development  and  structure  of  turbulent  flows  is  the  foous  of 
intense  study.  We  now  know  that  flows  which  previously  were  thought  to  be 
totally  chaotio  and  statistical  in  nature  are  dominated  by  the  persistence 
of  relatively  large  structures.  These  coherent  structures  were  ignored 
since  their  existence  was  masked  or  de-emphasized  by  experimental 
averaging  techniques.  The  classical  description  of  turbulence  and  the 
mechanisms  responsible  for  its  development  are  now  considered  deficient  in 
their  explanations  of  these  transient  but  organized  and  persistent 
structures  which  transport  much  of  the  stress  in  a  shear  layer. 

The  classical  description  of  turbulence  evolved  from  the  observed 
behavior  of  fluid  flows  as  a  function  of  Reynolds  number.  Many  flows 
exhibit  a  series  of  sudden  transitions  to  new  flow  patterns  as  the 
Reynolds  number  is  increased.  Each  transition  results  in  an  increasingly 
complicated  flow,  and  at  high  Reynolds  number  the  flows  become  irregular 
and  appear  chaotic  in  both  space  and  time.  The  transitions  to  the 
succession  of  flow  patterns  may  be  caused  by  a  sequence  of  fluid 
instabilities,  each  of  which  breaks  some  symmetry  in  the  previously  stable 
flow  pattern  and  introduces  some  new  scale  in  the  new  flow  pattern 
(Liepmann  1979).  The  transition  to  a  purely  chaotio  flow  was  postulated 
to  occur  through  an  infinite  succession  of  instabilities,  each 
contributing  to  the  increasing  frequency  content  of  the  flow  (Landau  & 
Lifshitz  1959). 

Several  experimental  observations  have  seriously  eroded  confidence  in 
the  completeness  of  such  a  turbulence  model  for  physioal  flows.  First, 
although  laboratory  created  grid  turbulenoe  comes  dose,  no  flow  has  yet 
been  shown  to  exhibit  pure,  homogeneous  isotropic  turbulence  in  the 
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classical  sense.  Such  a  state  is  really  a  limiting  condition.  Further, 
since  coherent  structures  exist  on  the  larger  scales,  it  is  reasonable  to 
assume  the  existence  of  similar  structures  on  all  scales  larger  than  the 
dissipation  scale.  Second,  the  transition  to  turbulence  does  not  occur 
through  an  infinite  succession  of  instabilities,  but  after  the  appearance 
of  relatively  few  Instabilities,  typically  three  or  four.  This  has  led  to 
the  postulate  of  a  theoretical  connection  between  strange  attractor  theory 
and  transition  to  turbulence.  Third,  intermittency  in  free  turbulent 
shear  layers  indicates  the  presence  of  a  thin,  sharp  interface  between 
turbulent  fluid  and  irrotational  fluid.  This  finding  leads  to  questions 
about  whether  such  sharp  interfaces  can  be  represented  as  a  diffusive 
effect.  Finally,  the  recent  dlsoovery  that  coherent  structures  dominate 
flows  which  were  previously  believed  to  approach  pure  turbulence  directly 
causes  us  to  reconsider  the  basic  assumptions  in  the  classical  theory  of 
turbulence. 

The  search  for  new  concepts  of  the  nature  of  turbulence  has 
centered  on  understanding  the  mechanisms  of  the  transition  to  turbulence 
in  several  simple  fluid  flows.  Shear  flows  generated  by  splitter-plate 
partitions  (e.g. ,  Brown  &  Roshko  1974,  Browand  &  Weidman  1976,  Roshko 
1976)  exhibit  most  of  the  troublesome  intricacies  associated  with 
the  transition  and  have  properties  that  are  relevant  to  practical 
applications.  Both  larger  and  smaller  scales  of  chaotic  fluid  motion 
develop  in  these  systems.  Flow  near  the  sharp  trailing  edge  of  the 
splitter  plate,  where  the  parallel  free-streams  first  meet,  may  initially 
be  laminar  for  even  high  Reynolds  numbers,  and  several  distinct 
transitions  are  easily  discerned  before  the  onset  of  an  apparently  fully 
turbulent  shear  layer.  It  is  still  quite  difficult,  however,  to  pinpoint 
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at  which  point  in  the  flow  the  label  'turbulent*  is  applicable.  Indeed, 
large  scale  coherent  structures  have  been  found  to  dominate  well 
downstream  in  flows  which  appear  chaotic  on  smaller  scales.  By  using 
similarity  arguments,  for  example,  we  see  that  the  splitter-plate  flow  is 
always  dominated  by  ever  larger  coherent  structures.  Through  flow 
visualization  experiments,  a  great  deal  of  insight  has  been  achieved  into 
the  exact  mechanisms  involved  in  the  individual  instabilities  as  well  as 
some  indication  of  the  sequence  of  appearance  of  various  scales  of  motion. 

The  shear  layer  generated  by  the  splitter  plate  is  unstable  to  the 
Kelvin-Helmholtz  instability.  Small  perturbations  in  the  flow  grow  into 
nonlinear  waves  which  break  and  roll  up,  transforming  the  original 
vorticity  of  the  shear  layer  into  Isolated  clumps.  The  primary  wavelength 
generated  by  the  instability  is  usually  that  of  the  fastest  growing  mode 
for  that  particular  geometry  or  of  some  Impressed  wavelength  determined  by 
boundary  conditions  or  initial  conditions.  Further  development  of  the 
shear  layer  proceeds  through  the  pairing  of  vortex  clumps,  a  process  which 
may  be  repeated  many  times  downstream.  One  of  the  effects  of  pairing  is 
to  generate  subharmonics  of  the  original  unstable  wavelength,  but  smaller 
wavelength  disturbances  are  oreated  as  well.  These  disparate  wavelengths 
arise  from  at  least  two  causes:  imperfect  pairing  due  to  small  fluctuations 
in  the  flow  leaving  an  unpaired  vortex  which  then  merges  with  a  previously 
formed  pair  (Browand  &  Winant  1973),  and  the  generation  of  small-scale 
disturbances  in  the  interaction  of  the  cores  of  the  two  vorticity  clumps 
(Zabusky  1981). 

A  number  of  experimental  investigations  have  shown  that  the  large 
spanwise  coherent  structures  dominate  the  entrainment  and  mixing  processes 
at  the  shear  layers  (Winant  A  Browand  1974,  Brown  &  Roshko  1974,  Dimotakis 


&  Broun  1976,  Breidenthal  1981).  Fluid  entrained  from  the  free  streams 
can  be  transported  unmixed  across  the  shear  layer.  When  chemical  reaction 
between  the  entrained  fluids  is  possible,  it  takes  place  preferentially  at 
the  regions  associated  with  the  coherent  structures.  In  many  practical 
applications,  such  as  combustion  and  chemical  lasers,  mixing  layers  are 
Important  and  it  is  extremely  helpful  to  understand  mechanisms  by  whloh 
mixing  can  be  enhanced  or  retarded. 

Recent  splitter-plate  experiments  by  Kooohesfahani  et  al.  (1983)  show 
that  more  high  speed  fluid  than  low  speed  fluid  is  entrained  in  the 
coherent  structures.  The  experimental  work  of  Konrad  (1976)  suggested 
that  the  entrainment  at  the  shear  layer  is  not  necessarily  symmetric. 
Broadwell  &  Breidenthal  (1982)  developed  a  simple  mixing  model  which  may 
explain  the  observed  results.  This  asymmetry  is  not  predicted  by  models 
based  on  gradient  transport  and  eddy  diffusivity.  It  has  also  not  been 
noted  in  spectral  simulations  of  unsteady  shear  flows  with  periodic 
boundary  conditions  (Riley  &  Metcalfe  1980),  or  in  previous  simulations  of 
the  splitter-plate  geometry  using  either  vortex  (Ashurst  1979)  or  finite 
difference  (Davis  &  Moore  1983)  methods. 

In  the  work  presented  below,  we  describe  numerical  simulations  of 


two-dimensional  shear  flows  of  a  splitter-plate  geometry  which  is  unstable 


to  the  Kelvin-Helmholtz  instability.  Our  goal  is  to  produce  accurate 


enough  calculations  of  the  time-evolution  of  the  major  physical  quantities 


so  that  we  can  analyze  the  flow  behavior  and  test  concepts  of  the 


transition  to  turbulence.  We  are  Interested  in  the  initial  transient 


behavior,  the  resulting  unsteady  pattern  of  structures,  and  the  symmetry 


properties  of  the  composition  of  their  associated  mixing  layers. 
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III.  THE  NUMERICAL  MODEL 


A.  Flux  Corrected  Transport 

The  FAST2D  simulation  model  was  used  to  perform  the  shear  flow 
calculations  described  below.  This  program  solves  the  compressible 
time-dependent  conservation  equations  for  mass,  momentum  and  energy 
coupled  to  algorithms  describing  gravity,  molecular  and  thermal  diffusion, 
and  chemical  reactions  with  energy  release.  The  different  algorithms  are 
coupled  by  time  step  splitting  so  the  corresponding  physical  effects  can 
be  turned  on  or  off  Independently  as  required  by  the  particular  problem 
being  studied. 

The  continuity  equations  are  solved  using  Flux-Corrected  Transport, 
(Boris  &  Book  1976),  a  non-linear,  monotone,  finite-difference  technique 
which  conserves  mass,  momentum  and  energy  and  preserves  the  physical 
positivity  property  of  fluids.  FCT  is  particularly  useful  in  problems 
where  sharp  discontinuities  arise  and  are  maintained  throughout  the 
calculations.  These  discontinuities  may  be  shocks  or  contact  surfaces, 
or,  as  in  the  case  presented  in  this  paper,  interfaces  ir  material 
density.  FCT  adds  linear  velocity-dependent  diffusion  to  a  high-order 
algorithm  during  convective  transport.  This  prevents  dispersive  ripples 
from  arising,  and  thus  ensures  that  all  conserved  quantities  remain 
monotone  and  positive.  Fourth  order  accuracy  is  maintained  by  subtracting 
out  this  added  diffusion  during  a  second  nonlinear  antidiffusion  stage  of 
the  algorithm  in  regions  where  it  is  not  needed  to  preserve  monotonicity. 
The  algorithm  has  been  used  extensively  in  one-,  two-  and  three-dimensional 
flows  (e.g.,  Oran  et  al.  (1982),  Picone  et  al.  (1984),  Kailasanath  et  al. 

( 1985a, b),  Fyfe  et  al.  (1985)). 


Another  important  feature  of  FCT  is  its  ability  to  separate  grid 
motions  from  the  fluid  flow.  This  allows  us  to  use  variably-spaced  grids 
as  well  as  moving  adaptive  grids  (Book  et  al.  1980,  Oran  et  al.  1982). 
Figure  1  is  an  example  of  a  convergence  test  that  was  done  for  the 
calculation  of  a  shock  which  transitioned  to  a  detonation.  A  shock 
propagating  in  a  tube  containing  a  mixture  of  hydrogen,  oxygen,  and  argon 
reflects  from  a  rigid  wall.  Upon  reflection,  the  temperature  and  pressure 
are  raised  high  enough  to  initiate  chemical  reactions.  The  chemical 
reactions  occur  first  near  the  wall,  where  a  reaction  wave  eventually 
forms.  The  reaction  wave  overtakes  the  reflected  shock  and  the  result  is 
a  propagating  detonation.  '  This  problem  was  studied  computationally  using 
FCT  coupled  to  a  chemical  kinetics  model,  and  the  results  were  compared  to 
experimental  shock  tube  data  (Oran  et  al.,  1982).  Figure  1  shows  a  series 
of  resolution  tests  for  one  of  the  cases  studied.  The  first  calculation, 
(a),  used  a  uniform,  stationary  grid  with  a  cell  spacing  of  0.07  cm.  The 
second  and  third  calculations,  (b)  and  (c),  both  used  moving,  variably 
spaced  grids  with  different  minimum  spacing  around  the  detonation  front. 
This  was  a  difficult  problem  to  formulate  because  the  grid  was  both 
variable  and  moving.  However,  the  cost-benefits  of  using  this  flexible 
grids  are  substantial.  Various  resolutions  tests  for  multidimensional 
versions  of  this  model  coupled  to  variably  spaced  and  moving  grids,  have 
been  reported  by,  e.g.,  Book  et  al.  (1980)  and  Kailasanath  et  al. 

( 1985a, b). 

The  particular  version  of  FCT  used  here  is  JPBFCT,  which  is  similar 
to  ETBFCT  (Boris  1976a).  This  routine  solves  the  one-dimensional 
continuity  equation  in  Cartesian,  cylindrical,  spherical  or  generalized 
nozzle  geometries,  depending  on  the  value  of  a  logical  switch.  Since  the 


algorithm  is  one-dimensional,  time  step  splitting  in  the  various 
directions  is  used  to  construct  two-  and  three-dimensional  codes.  The 


two-dimensional  Cartesian  calculations  described  below  were  performed  with 
the  general  code  with  the  gravity,  diffusion,  chemistry,  and  energy 
release  options  not  used  (discussed  elsewhere  by  Oran  4  Boris  1981). 

B.  Computational  Grid 

The  computational  grid  was  set  up  initially  in  the  calculations  and 
then  held  fixed  in  time.  A  typical  grid  is  shown  in  figure  2(a).  Finely 
spaced  cells  clustered  around  the  centerline  where  the  instability  first 
occurs  and  the  coherent  structures  form.  The  splitter  plate  is  on  the 
centerline  and  its  trailing  edge  is  at  T  =  2.0  cm.  The  Y-direction  is 
well  resolved  near  the  centerline.  The  velocity  of  the  lower  stream  (Vf) 
is  always  considerably  larger  than  that  of  the  upper  stream  (V  ).  As  the 
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large  spanwise  vortices  propagate  downstream,  they  merge  and  grow.  The 
idea  is  to  have  the  larger  computational  cells  further  downstream  and 
further  from  the  centerline  resolve  the  larger  coherent  structures,  and 
the  smaller  cells  resolve  the  smaller  structures. 

Grid  spacing  was  varied  in  different  calculations  to  (1)  check  the 
effects  of  resolution,  (2)  test  the  implementations  of  various  boundary 
condition  algorithms,  and  (3)  test  the  effects  of  physically  different 
boundaries.  Maximums  of  100  cells  across  the  mixing  layer  (X-direction) 
and  170  cells  in  the  streamwise  direction  (Y-direction)  were  used.  The 
first  20  of  the  cells  in  the  Y-direction  were  in  the  region  behind  the  tip 
of  the  splitter  plate.  The  mesh  spacing  varied  in  the  range  0.024  -  0.24 
cm  in  the  X-direction  and  from  0.06  -  0.46  cm  in  the  Y-direction.  In 
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general,  the  spacing  should  not  stretch  or  compress  more  than  20-30%  from 
cell  to  cell. 

Typical  convergence  tests  varied  the  spacing  in  both  directions  by 
plus  or  minus  50$.  The  observed  differences  in  the  calculations  made 
minor  differences  in  the  quantitative  results,  and  no  differences  in  the 
qualitative  observations.  In  an  Eulerian  finite-difference  calculation, 
varying  the  grid  spacing  is  the  most  rigorous  test  of  a  computational 
model.  The  results  presented  in  Section  IV  have  been  calculated  on  a 
60x150  grid.  A  discussion  of  the  sensitivity  of  the  results  to  the 
location  of  the  boundaries  and  the  boundary  conditions  used  in  the 
calculations  is  given  at  the  end  of  Section  IV. 

C.  Boundary  Conditions 

As  shown  in  figure  2(b),  the  top  and  the  bottom  of  the  computational 
region  are  perfectly  reflecting  hard  walls.  These  boundary  conditions  are 
implemented  by  assuming  that  the  X-component  of  the  velocity  at  the  guard 
cell  is  equal  to  minus  the  X-component  of  the  velocity  at  the  first  cell 
inside  the  domain.  These  type  of  boundary  conditions  are  also  imposed  at 
the  splitter  plate  to  ensure  that  there  is  no  flux  of  material  aoross  its 
surface.  Assuming  a  free-slip  condition  at  the  walls  is  a  first  approximation 
that  we  use  in  this  inviscid  calculation.  It  is  unlikely  that  the  effects 
of  boundary  layers  at  the  walls  would  considerably  affect  the  calculated 
flowfield.  In  particular,  the  location  and  presence  of  the  walls  does  not 
significantly  affeot  the  calculated  asymmetries  described  below  in  Section  III.C 

The  fundamental  difficulty  when  solving  problems  with  outflow 
boundary  conditions  is  that  information  about  the  flow  beyond  the 
computational  mesh  is  required  to  make  the  fluid  near  the  boundaries 


behave  properly.  This  problem  is  generally  handled  by  using  guard  cells. 
Guard  cells  are  just  outside  the  computational  regime,  and  are  not 
actually  part  of  the  calculation.  They  are  used  to  tell  the  boundary 
cells  how  the  outside  world  is  behaving.  The  simplest  model  of  outflow  in 
guard  cells  is  to  say  that  the  momentum,  energy,  and  density  do  not 
change,  i.e.,  there  is  effectively  zero  gradient.  This  causes  problems  in 
extended  calculations  since  it  does  not  provide  a  way  for  the  solutions  to 
relax  to  background  conditions. 

An  outflow  boundary  algorithm  was  developed  to  use  with  the  FCT 
algorithm  described  in  Section  III. A.  This  algorithm  defines  the  values  of 
quantities  at  the  guard  cell  from  a  zeroth-order  extrapolation  of  the 
value  at  the  boundary  cell.  In  addition,  it  assumes  that  there  is  a  slow 
local  relaxation  outside  the  mesh  towards  the  known  ambient  value.  The 
strong  nonlinear  stabilizing  properties  of  the  FCT  method  appear  to 
eliminate  instabilities  which  occur  in  other  nonlocal  methods  when  low 
order  extrapolations  are  used  for  specifying  boundary  conditions  (Turkel, 
1980).  Previous  tests  (Boris  et  al.  1985)  have  suggested  this  as  a 
simple,  yet  effective  approach  to  the  outflow  boundary  conditions.  Such 
an  expression  is  expected  to  be  an  approximation  to  the  lowest  order  terms 
in  an  asymptotic  expansion,  valid  for  times  short  compared  to  a  sonic 
transit  time  of  the  system. 

For  two-dimensional  incompressible  calculations,  Davis  &  Moore  (1983) 
specify  background  values  and  use  zero-gradient  conditions  at  infinite 
distances  from  the  computational  region.  Their  boundary  conditions  are 
Implemented  by  mapping  the  infinite  physical  space  into  a  finite 
computational  domain.  Their  approach  to  outflow  has  important  elements  in 
common  with  the  method  used  here.  Both  methods  stretoh  the  mesh  to  get 


the  boundary  far  from  the  region  of  interest. 

The  algorithm  is  implemented  by  assuming  that  the  pressure  in  the 
guard  cell  is  extrapolated  from  the  boundary  cell  but  relaxes  to  a 
prespecified  ambient  value.  Specifically,  if  7  is  the  streamirlse 
(outflow)  direction. 


(7  -7  ) 

P  =  Pn  +  “  P„> 

g  n  (y  ▼  )  aBb  n 

g"  s' 
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where  is  the  guard  cell  pressure,  the  ambient  pressure^  Pq  is  the 

boundary  pressure  and  7  ,  7  and  7  ,  correspond  to  the  positions  of  the  guard 
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cell,  the  boundary  cell  and  the  cell  at  the  edge  of  the  splitter  plate, 
respectively.  This  expression  is  derived  by  assuming  that  the  pressure 
reaches  its  ambient  value  at  7  *  co  ,  and  P  is  then  obtained  by  linearly 

o 

Interpolating  (in  the  variable  1/(7-7  )  )  between  P  and  P  .  .  This  is  a  one 

s  n  amb 

dimensional  first  approximation  to  the  correct  outflow  boundary  condition. 

The  density  and  momenta  are  extrapolated  according  to 
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Further  improvements  for  Pg  are  possible  by  considering  higher  order 
Interpolation  schemes,  i.e.,  by  including  the  information  on  the  pressure 
at  other  cells  in  the  domain  in  addition  to  that  at  the  boundary  cells. 

In  particular,  it  might  be  useful  to  include  an  explicit  dependenoe  on 
physical  parameters  of  the  problem  such  as  the  separation  between  the 
walls,  the  size  and  rate  of  flow  of  the  structures  and  the  local  speed  of 
sound.  This  requires  further  physical  assumptions  about  conditions 
outside  of  the  computational  domain  in  the  downstream  direction. 
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In  our  preliminary  calculations,  ve  used  Inflow  boundary  conditions 
wtalch  prespecified,  constant  values  of  the  mass,  momentum  and  energy  of 
the  Inflowing  gas.  Also,  the  inflow  boundary  was  at  the  tip  of  the 
splitter  plate.  These  conditions,  however,  did  cot  correctly  provide  the 
feedback  between  the  fluid  just  entering  the  computational  domain  and  the 
disturbances  created  downstream  by  the  Kelvin-Helmholtz  roll-ups  and 
vortex  merging.  Pressure  pulses  from  downstream  disturbances  create  small 
transverse  flows  at  the  trailing  edge  of  the  splitter  plate.  These  pulses 
re-initiate  the  instability  and  lead  to  the  next  ooherent  vortex  roll  up. 
When  the  inflow  pressure  was  held  equal  to  the  ambient  pressure,  the  first 
vortex  structures  formed  very  far  downstream.  They  eventually  flowed  off 
the  computational  domain,  and  the  instability  appeared  to  die  out.  The 
apparent  reason  for  this  is  that  the  pressure  perturbations  arriving  at 
the  inflow  boundary  were  cancelled  by  the  non-physical  condition  whioh 
kept  the  inflow  pressure  constant.  A  more  physically  reasonable  treatment 
of  the  Inflow  response  to  pressure  fluctuations  was  required. 

The  Inflow  boundary  conditions  we  use  now  specify  the  inflow  density 
and  velocity,  and  then  use  a  zero-slope  condition  on  the  pressure  at  the 
inflow  boundary  to  derive  the  energy.  This  algorithm  allows  pressure 
differences  between  the  top  and  bottom  streams  to  generate  transverse 
flows.  In  addition,  a  short  inflow  plenum  is  modelled  by  including  the 
2.0  cm  region  behind  the  end  of  the  splitter  plate,  in  the  computational 
domain,  as  shown  in  figure  2(a). 

With  a  fluctuating  Inflow  pressure  condition,  it  is  important  to 
relax  the  outflow  pressure  toward  an  ambient  value.  This  is  because  a 
base  pressure  for  the  system  has  to  be  specified  in  compressible 
calculations  when  the  value  of  the  Inflow  pressure  is  allowed  to  vary. 


Relaxation  of  the  density  and  momentum  towards  their  ambient  values  is 
less  Important  on  outflow  because  these  values  are  given  at  the  inflow 
boundary.  The  long  time  values  of  all  the  variables  must  be  available  to 
the  solution  from  the  beginning  to  prevent  secular  deviation  of  the 
calculated  flow. 


IV.  SPLITTER-PLATE  SIMULATIONS 


As  shown  in  figure  2(b),  two  coflowing  streams  of  air  are  Initially 
separated  by  a  thin  plate  and  then  enter  a  long  chamber.  The  simulations 
are  Initialized  by  assuming  that  both  gas  streams  have  the  same  initial 

pressure  (1  atm)  and  temperature  (298  K),  and  have  initial  uniform 

velocities  V.  and  V  .  The  values  considered  for  V.  and  V  ranged  from 

1  3  1  8 

b  b  A 

1.0  z  104  to  2.0  x  10*  cm/s,  and  from  0.0  to  4.0  x  ICr  cm/s,  respectively. 

The  range  V  /V.  of  0.0  -  0.2  was  studied.  Sinoe  the  velocity  of  sound  in 

s  r 
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this  mixture  is  typically  of  the  order  of  3.5  x  10  cm/s,  the  flows  were 
subsonic  and  hence  virtually  incompressible.  The  PCT  algorithm  used  to 
describe  the  convection  flow,  however,  is  fully  compressible  so  sound 
waves  are  resolved  in  the  system  and  acoustic  delay  times  for  pressure 
waves  are  properly  accounted  for. 

The  small  perturbation  which  initiates  the  transition  from  laminar  to 
turbulent  flow  oocurs  at  the  first  tlmestep  of  the  calculation.  It  is 
analogous  to  using  a  delta  function  perturbation  at  the  center  in  a 
perlodlo  simulation  of  two  equal  and  opposite  streams.  The  perturbation 
generates  small  pressure  gradients  and  diffuses  vorticity  at  the  shear 
layer,  near  the  edge  of  the  splitter  plate.  The  initial  step  funotlon 
profile  for  the  downstream  velocity  is  smootbed-out  after  the  first  few 
convection  cycles.  As  a  consequence,  the  initially  uniform  pressure 
develops  a  small  gradient  across  the  tbiokening  shear  layer  in  the  region 
close  to  the  edge  of  the  splitter  plate  where  the  free  streams  meet.  This 
disturbance  moves  downstream  as  the  integration  proceeds,  generating  the 
transverse  flows  which  trigger  the  instability. 
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A.  Transition  from  Laminar  to  Turbulent  Flow. 
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Figures  3-7  show  the  development  of  the  flow  for  the  case  in  which 

V,  =  1.0  x  104  cm/s  and  V.  =  2.0  x  103  cm/s  (i.e.,  V  /V_  =  0.2  ). 
r  s  3  r 

The  time  evolution  of  the  flow  is  shown  through  sequences  of  isovorticity 
and  number  density  ratio  (R)  contour  plots,  where 


R  =  . 

1TTF 


(5) 


The  and  N8  are  the  number  densities  corresponding  to  the  faster  and 
slower  inflowing  streams,  respectively.  R  varies  between  0  and  1,  which 
correspond  to  100*  of  the  slower  stream  or  10 Of  of  the  faster  stream, 
respectively. 

In  figures  3  and  5,  H  is  contoured  in  the  interval  0.3  -  0.7.  This 
particular  choice  of  contour  levels  is  useful  for  describing  the  mixing. 
The  transition  from  a  uniform  shear  flow  is  first  noticeable  by  the 
appeareance  of  a  pair  of  vortioes  forming  just  ahead  of  the  tip  of  the 
plate,  shown  at  0.525  ms  In  figure  3.  As  this  pair  travels  downstream,  it 
grows  steadily  due  to  the  vorticity  whloh  is  being  fed  in  to  it  from  both 
directions.  In  addition,  smaller  structures  are  entrained  by  the  roll-up 
on  the  side  of  the  splitter-plate.  As  the  system  evolves,  the  smaller 
struotures  have  more  time  to  grow  before  they  are  engulfed.  The  last 
panel,  at  0.975  ms,  shows  one  of  these  small  structures,  now  a  vortex, 
being  shed  and  merging  with  the  large  vortex. 

The  isovorticity  contours  in  figure  4  show  the  shear  layer  breaking 
up.  New  vorticity  dumps  form  and  subsequently  pair  with  the  leading 
vortex.  A  relation  exists  between  the  growth  of  the  leading  vortex, 


16 


through  the  pairing  process,  and  the  rate  at  which  new  vortices  are  shed. 
Pressure  disturbances  are  generated  by  the  fluid  accelerations  caused  by 
downstream  roll-ups  and  mergings.  These  disturbances  can  reinitiate  the 
instabilities  close  to  the  trailing  edge  of  the  splitter  plate.  The 
choice  of  the  inflow  boundary  condition  was  crucial  in  allowing  for  this 
physical  feedback  mechanism  to  occur.  This  concept  was  tested  by  damping 
the  acoustic  pulses  reaching  the  inflow  boundary.  This  was  the  Initial 
inflow  boundary  condition  described  above  in  Section  III,  where  the 
pressure  at  the  tip  of  the  splitter  plate  was  set  equal  to  the  ambient 
pressure. 

Figure  5  is  a  continuation  of  the  same  calculation  at  later  times, 
from  0.975  to  2.620  ms.  Small  vortices  are  shed  at  about  1.0  -  1.6  cm 
from  the  tip  of  the  splitter  plate.  Larger  vortices  have  moved 
downstream.  Pairings  between  the  newly  shed  vortices  are  also  possible 
now,  as  they  become  more  separated  from  the  leading  vortex.  By  2.62  ms, 
the  last  panel  in  figure  5,  merging  vortices  move  off  the  right  edge  of 
the  computational  system  and  unperturbed  fluid  moves  in  on  the  left. 

Figure  6  shows  contours  of  R  in  the  range  0.70  -  0.97.  This  choice 
of  contour  levels  emphasizes  the  portions  of  the  faster  fluid  which  are 
transported  largely  unmixed  across  the  shear  layer.  Isovorticity  contours 
corresponding  to  these  frames  are  shown  in  figure  7.  We  see  that  the 
growth  pattern  of  the  structures  appears  to  follow  a  similarity  law.  This 
is  consistent  with  the  flow  having  a  length  scale  which  increases 
downstream  proportionally  to  the  distance  from  an  origin  of  similarity 
located  just  ahead  of  the  edge  of  the  splitter  plate.  By  the  last  frame 
at  2.62  ms  (see  also  figure  5),  a  spectrum  of  modes  is  seen  suggesting  a 
snapshot  from  an  experimental  flow  visualization.  As  can  be  seen  from  the 


figures,  some  asymmetry  exists  In  the  similarity  pattern.  Structures 
penetrate  further  from  the  center  line  on  the  low  speed  side.  This 
deflection  is  also  observed  in  the  experiments  (see,  e.g.,  the 
shadowgraphs  in  Brown  &  Roshko  (1974)  and  in  Konrad  (1976)). 

B.  Asymmetries  in  the  Composition  of  the  Mixing  Layers. 

Experiments  by  Koochesfahani  et  al.  (1983)  show  that  there  is  more 
high  speed  than  low  speed  material  entrained  and  mixed  in  the  coherent 
structures.  These  experiments  investigated  the  composition  of  the  mixing 
layers  by  using  a  dilute  acid  solution  and  a  dilute  base  solution  for 
each  of  the  free  streams,  and  examining  the  extent  of  the  titration 
process.  A  pH  sensitive  fluorescent  dye  was  premixed  with  the  acid 
solution.  Then  laser  induced  fluorescence  was  used  to  monitor  the  local 
chemical  environment  at  the  shear  layer  and  in  the  mixed  region  between 
tne  two  streams.  The  amount  of  chemical  product  was  estimated  as  a 
function  of  the  amount  of  (fluorescent)  fluid  having  a  pH  above  the 
fluorescence  threshold.  Two  Identical  experiments  were  performed  in  which 
the  faster  stream  was  alternately  acidic  or  basic.  Comparison  of  the 
results  indicated  that  the  amount  of  product  was  larger  by  two  orders  of 
magnitude  in  the  case  in  whioh  the  faster  stream  was  basic..  This  suggests 
that  the  entrainment  at  the  shear  layer  is  asymmetric,  with  an  excess  of 
higher  speed  fluid  in  the  composition  of  the  mixing  layers. 

In  order  to  evaluate  the  mixing  asymmetries  in  our  calculations,  we 
have  simulated  this  fluorescence  diagnostic.  We  assume  that  the  streams 
of  air  are  prepared  as  in  the  experiment,  namely,  that  one  is  a  dilute 
'acid*  solution  carrying  a  pH  sensitive  dye  whioh  ’fluoresces'  when  the  pH 
is  above  a  certain  threshold,  and  the  other  is  'basic'.  The  add  (plus 


dye)  and  base  additives  to  the  streams  of  gas  are  small  enough  that 
eonveotion  is  not  perturbed  by  their  presence. 


We  then  define  a  ’ fluorescence ’  intensity,  I.  Let  streams  's'  and 
'f'  be  the  'acid'  and  'base'  streams,  respectively.  Also  let  the  streams 
have  an  equivalence  ratio  such  that  1  part  of  the  'base'  stream  is 


required  to  neutralize  1  part  of  the  'acid'  stream  and  elevate  the  local 
pH  above  the  dye  threshold.  Then,  the  'fluorescence'  intensity  at 
at  cell  (i,j)  is 


(6) 


where  N  is  a  normalization  constant.  The  index  'i'  refers  to 
o 

the  cross-stream  (X)  direction  and  '  J '  refers  to  the  streamwise  (T) 
direction.  In  addition,  we  assume  that  a  'fluorescent'  portion  of  the 
fluid  has  an  excess  of  base,  as  in  the  experiment,  and  that  the  intensity 
is  proportional  to  the  amount  of  dye  in  the  cell.  This  ensures  that  the 
Intensity  decreases  to  zero  in  the  pure  base  stream  and  provides  a  very 
simple  physical  cut-off  defining  the  mixing  layers.  Analogously,  we 
define  an  intensity  I8^^  corresponding  to  the  reversed  case  in  which 
streams  'f'  and  's'  correspond  to  the  'add  +  dye'  and  'base*  solutions, 
respectively.  Both  intensities  can  then  be  calculated  from  the  same 


computed  data. 

By  comparing  I*"8  and  Isf,  we  can  evaluate  the  mixing  asymmetry  at 
the  shear  layer  in  terms  similar  to  those  used  in  the  experiment.  A 
smooth  exponential  interpolation  was  made  between  0.0  and  in 

defining  for  in  the  range  0.95  x  N8^  to  1.0  x  H8^.  The 

corresponding  Interpolation  was  made  for  I8**^ .  The  normalization 


constant  Nq  was  defined  as  Nq  =  P/k^T,  where  P  and  T  are  the  values  of 

the  initial  uniform  pressure  and  tempature,  and  kg  is  Boltzmann's 

constant.  With  this  definition,  the  'fluorescence'  intensities  are 

between  0.0  and  0.5.  Some  representative  cases  are  shown  in  figures 

8-10  for  7  /V  s  0.0,  0.1,  and  0.2.  The  first  case  corresponds  to  a 
s  r 

quiescent  upper  stream.  The  colored  regions  indicate  the  locations 
where  the  intensities  are  appreciably  different  from  zero,  as  determined 
by  the  requirement  I  >  0.005. 

Figures  8-10  show  that  as  the  velocity  ratio  increases,  the 
instability  appears  closer  to  the  edge  of  the  splitter  plate,  and  the 
structures  grow  faster.  This  result  is  expected,  since  the  growth  rates 
of  the  Kelvin-Helmholtz  instability  increase  with  the  velocity  difference. 
We  can  roughly  define  fluorescent  regions  as  those  regions  covered  by  the 
lighter  colors,  i.e.,  those  for  which  I  is  greater  than  a  threshold  value 
Iq  s  0.125,  equal  to  25$  of  the  maximum  attainable  value.  The  figures  show 
that  the  fluorescent  area  is  greater  when  the  faster  stream  is  basic  and 
the  dye  is  in  the  slower  acid  stream.  This  means  that  the  mixing  layer 
has  more  of  the  faster  than  of  the  slower  fluid. 

In  order  to  gat  a  global  measure  of  the  asymmetry,  we  consider  a 
ratio  of  effective  fluorescent  areas,  Q(Y),  defined  by 
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The  sum  over  the  index  'i'  is  in  the  cross-stream  X-direction  and  is 


the  number  of  cells  in  that  direction.  The  sum  over  the  index  ' J'  extends 


over  computational  cells  from  the  tip  of  the  splitter  plate,  at  location 

T  ,  cell  J  ,  to  a  cut-off  distance  7  at  cell  J.  The  quantity  A. .  is  the 
o'  o  ij 

area  of  the  cell  (i,j). 

The  numerator  and  denominator  in  equation  (7)  correspond , respectively 
to  the  fluorescent  areas  for  the  "fs"  (base  in  the  faster  fluid,  acid  + 
dye  in  the  slower  fluid)  and  ”sf”  (acid  +  dye  in  the  faster  fluid  and  base 
in  the  slower  fluid)  flow  configurations.  Then,  a  mixing  region  of  length 
(T  -  Yq)  contains  more  of  the  faster  fluid  if  Q(T)  is  greater  than  1.0. 

The  average  magnitude  of  (Q  -  1.0)  in  the  interval  (Yq,  7)  considered 
gives  a  measure  of  the  mixing  asymmetry.  Figure  11  shows  Q(Y)  vs.  7  for 

4 

two  cases,  V  /V.  =0.1,  and  0.2,  with  V-  *  2  x  10  cm/s  at  0.9  ms.  This 
9  r  r 

figure  also  shows  that  the  mixing  asymmetry  increases  as  the  ratio  of 
free-streem  velocities  Increases. 

C.  The  Effects  of  Boundaries  on  the  Mixing  Asymmetry 

We  performed  a  number  of  tests  to  evaluate  the  sensitivity  of  the 
calculations  to  the  location  of  the  boundaries  and  to  different  forms  of 
outflow  boundary  conditions.  Our  studies  Indicated  that  the  results  were 

i 

not  particularly  sensitive  to  these  kinds  of  variations.  Some 
representative  results  involving  rather  significant  changes  are  shown 
below  in  figures  12  and  13. 

To  test  the  effects  of  the  location  of  the  boundaries,  we  added 

computational  cells  above  and  below  the  splitter-plate  in  the  cross-stream 

direction,  and  at  the  outflow  boundary  in  the  streamwise  direction. 

Figure  12  shows  the  ratio  of  fluorescent  areas,  Q,  for  two  different  sized 

4 

computational  regions  for  the  case  in  which  Vs/Vf  =  0.1,  Vf  =  2.0  x  10 

p 

cm/s,  at  0.960  ms.  The  smaller  region  is  the  6.13  x  19.44  cm  considered 

2 

before,  and  the  other  is  a  larger  22,30  x  24.9  cm  domain.  The  shapes  and 


exact  locations  of  the  maxima  and  minima  of  the  curves  are  somewhat 
different  since  the  curves  are  not  for  physically  identical  problems. 
However,  the  general  features  and,  in  particular,  the  asymmetry  patterns, 
are  essentially  the  same. 

In  figure  13  we  show  the  sensitivity  of  the  results  to  changes  in  the 
outflow  boundary  condition.  We  compare  the  curves  of  Q  for  two  different 
outflow  boundary  conditions  in  the  case  where  V  /V.  =  0.2,  V_  =  1.0  x  10** 
cm/s,  at  1.725  ms.  The  solid  line  is  for  an  outflow  boundary  condition 
in  which  the  pressure  at  the  outflow  boundary  relaxes  to  its  background 
value,  as  given  in  equation  (1).  The  solid  line  with  circles  is  for 
calculations  in  which  the  streamwise  pressure  gradient  was  set  equal  to 
zero  at  the  outflow.  In  this  case  (zero-th  order  extrapolation)  we 
neglected  the  second  term  on  the  right  hand  side  of  equation  ( 1 ) . 

Although  this  condition  is  unphysical,  it  gives  a  useful  limiting  case. 

For  the  times  considered,  the  differences  between  the  curves  are  only 
slight,  showing  that  the  results  are  insensitive  to  changes  in  the 
outflow  boundary  conditions.  The  pattern  of  asymmetry  is  essentially 
the  same.  The  figure  shows,  however,  that  the  effect  of  the  zero-th 
order  pressure  extrapolation  condition  at  this  stage  is  to  tend  to 
smooth  the  structures. 
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V.  DISCUSSION  AND  SUMMARY 


This  paper  has  presented  numerical  simulations  of  the  evolution  of 
the  Kelvin-Helmholtz  instability  for  a  splitter-plate  geometry.  These 
calculations  were  performed  using  the  Eulerian,  explicit  Flux-Corrected 
Transport  algorithm  to  solve  the  convective  transport  problem.  New 
boundary  conditions  were  developed  to  allow  the  inflowing  material  to 
respond  to  pressure  fluctuations  occurring  in  the  computational  region. 

New  boundary  conditions  were  also  developed  for  the  outflow  to  model  the 
relaxation  of  the  gas  in  the  computational  domain  to  the  unperturbed 
ambient  gas. 

The  focus  of  this  paper  is  the  evolution  of  the  coherent  structures 
in  spatially-evolving  mixing  layers,  with  particular  emphasis  on  their 
composition.  The  computations  show  a  trend  in  which  more  fast  fluid  than 
slow  fluid  is  entrained  in  the  st rue t\ires.  This  asymmetry  is  larger  as 
the  velocity  ratio  of  the  two  streams  increases.  The  calculated  effect 
is  consistent  with  the  experimental  results  of  Koochesfahani  et  al.  (1983), 
discussed  above.  We  note,  however,  that  the  experiments  were  done  using 
liquids,  while  the  calculations  considered  a  fast-moving  (but  subsonic) 
gas.  In  addition,  the  computational  model  only  measures  the  large  scale 
features  of  the  mixing  layers,  i.e.,  convective  mixing,  since  molecular 
diffusion  is  not  included. 

The  mixing  asymmetry  in  spatially  evolving  mixing  layers  arises  because 
of  a  combination  of  conditions.  First,  there  is  the  velocity  difference  and 
the  direction  of  circulation.  Added  to  this,  the  larger,  downstream  structures 
draw  in  more  faster  than  slower  fluid.  This  combination  of  effects  allows 
the  asymmetry  to  develop.  The  mixing  asymmetry  cannot  occur  in  simulations 
with  periodic  boundary  conditions  and  equal  and  opposite  velocities. 


The  co-flowing  system  of  two  streams  of  air  described  above  is 
unstable  as  long  as  the  Reynolds  number  is  high  enough  and  there  is  some 
small  perturbation  in  the  system.  This  perturbation  can  be  an  imposed 
frequency  at  the  beginning  of  the  calculation,  or  simply  noise  or  roundoff 
error.  Because  feedback  is  allowed  between  the  downstream  events  and  the 
inflow  boundary  conditions  through  the  zero  pressure  gradient  condition 
at  the  inflow,  there  is  no  need  to  constantly  drive  the  instability. 

The  mechanism  which  reinitiates  the  instabilities  close  to  the  inflow 
boundary  works  through  pressure  pulses  generated  at  various  scales  by  the 
fluid  accelerations  caused  by  downstream  vortex  rollups  and  mergings. 

They  are  required  to  ensure  that  the  flow  remains  essentially  divergence 
free  everywhere  when  a  coherent  structure  or  vortex  is  locally 
accelerated.  Though  these  pulses  are  transmitted  acoustically,  they  exist 
even  in  the  incompressible  limit. 

Evidence  supporting  this  mechanism  is  given  by  the  sequences  of 
panels  shown  in  figures  4  and  5.  The  flow  structure  at  0.975  ms,  for 
example,  is  both  large  and  irregular  enough  to  generate  pressure 
perturbations,  which  drive  noticeable  fluctuations  throughout  the 
computational  region.  In  addition,  the  reflection  of  the  pressure  pulses 
on  the  confining  walls  affects  the  development  of  the  larger  structures 
which,  in  turn,  influence  the  pattern  of  formation  of  new  structures 
upstream.  A  feedback  mechanism,  in  which  the  upstream  flow  is  influenced 
by  the  downstream,  was  proposed  by  Dimotakis  &  Brown  (1976).  This  was 
used  to  explain  unusually  long  autocorrelation  times  of  the  streamwlse 
velocity  fluctuations  observed  in  the  experiments.  These  could  not  be 
explained  by  a  local  flow  property. 
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A  number  of  tests  were  made  to  evaluate  the  Importance  of  the 
resolution  and  size  of  the  computational  region  and  the  location  of  the 
solid  bounding  walls  to  the  mixing  asymmetry.  The  resolution  tests  showed 
that  the  large-scale  features  and  the  mixing  asymmetry  reported  above  were 
adequately  resolved  by  the  computational  grid  chosen.  The  locations  of 
the  bounding  walls  did  not  appear  to  affect  the  basic  calculated  mixing 
asymmetry.  Also,  the  mixing  asymmetry  results  were  not  particularly 
sensitive  to  the  choice  of  the  outflow  boundary  conditions,  although  using 
a  physically  correct  inflow  boundary  condition  is  extremely  important. 

The  accuracy  of  the  FCT  calculations  on  variably-spaced  grids  has 
been  tested  extensively.  This  was  discussed  in  Section  III. A.  It  is 
important,  however,  to  consider  what  the  variable  grid  might  do  for  the 
splitter-plate  calculations  presented  here.  We  know  that  as  the 
structures  in  the  calculations  move  into  more  widely  spaced  zones, 
short-wavelength  information  is  damped.  Also,  as  the  spectrum  of  pressure 
disturbances  created  by  roll-ups  and  mergings  moves  into  more  finely 
resolved  regions  from  coarser  regions,  it  is  missing  a  portion  of  its 
short  wavelength  spectrum  that  it  would  have  if  the  zones  downstream  were 
finer.  The  result  we  have  seen  is  that  the  large-scale  structures  and 
features  of  the  asymmetry  are  not  significantly  affected  by  the  grid 
spacings.  Future  investigations  will  look  at  the  finer  details  of  how  the 
re-initlatlon  process  might  be  affected  by  these  missing  short  wavelengths 
propagating  from  downstream. 

The  exact  structure  of  the  unstable  flow  at  a  given  time  is  very 
sensitive  to  background  conditions,  initial  conditions,  and  to  small 
fluctuations  in  the  system.  This  is  the  crux  of  the  problem  in 
experiments  also:  many  realizations  are  possible  depending  on  fluctuations 
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in  the  initial  and  boundary  conditions.  Thus,  we  observe  that  some  of  the 
difficulty  encountered  in  the  computations  reflects  the  same  difficulty  in 
the  ezperi-  ments.  Both  highlight  the  highly  nonlinear  sensitivity  of  the 
system  to  the  boundary  and  initial  conditions. 
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Contours  of  R,  as  in  figure  3,  for  later  stages  in  the  development 
of  the  flow. 
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Contours  of  the  nixing  intensity  I,  as  defined  by  equation  6 
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for  V  /V-  =  0.0,  with  V-  =  2.0  x  10"  cn/s. 
s  r  i 

Dark  blue:  0.005  <  I  <  0.125,  Turquoise:  0.125  <  I  <  0.250, 
Green:  0.250  <  I  <  0.375,  Yellow:  0.375  <  I  <  0.500. 
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Figure  11.  Ratio  of  fluorescent  areas  as  a  function  of  the  streamwise 

a 

direction  at  0.9  ms,  for  Vf  =  2.0  x  10  cm/s; 
solid  line  :  Vs/Vf  =0.1,  solid  line  with  circles  : 


Ratio  of  fluorescent  areas  as  a  function  of  the  streamwise 

direction  at  0.96  ms,  for  7  /V-  *  0.1  and  V-  *  2.0  z  101* 

9  r  r 

2 

cm/s;  solid  line  :  6.13  x  19.44  cm  domain,  solid  line  with 

2 

circles:  22.3  x  24.9  om  domain. 


Figure  13*  Ratio  of  fluorescent  areas  as  a  function  of  the  streamwise 

a 

direction  at  1.73  ns,  for  V  /V,  =  0.2  and  V.  =  1.0  x  10 

Si  i 

cm/s;  solid  line  :  relaxation  of  outflow  boundary  pressure 
using  equation  (1),  solid  line  with  circles  :  zero-th  order 
extrapolation  of  the  outflow  boundary  pressure. 
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